Superfluidity without symmetry-breaking: the time-dependent Hartree-Fock 
approximation for Bose-condensed Systems 
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We develop a time-dependent Hartree-Fock approximation that is appropriate for Bose-condensed 
systems. Defining a depletion Green's function allows the construction of condensate and depletion 
particle densities from eigenstates of a single time-dependent Hamiltonian, guaranteeing that our 
approach is a conserving approximation. The poles of this Green's function yield the energies of 
number-changing excitations for which the condensate particle number is held fixed, which we show 
has a gapped spectrum in the superfluid state. The linearized time-dependent version of this has 
poles at the collective frequencies of the system, yielding the expected zero sound mode for a uniform 
infinite system. We show how the approximations may be expressed in a general linear response 
formalism. 
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I. INTRODUCTION 

The time-dependent Hartree-Fock approximation 
(TDHFA) is a powerful method for determining the prop- 
erties of many-body systems 0. It allows the computa- 
tion of the density response functions, conductivities, and 
other linear response functions. The poles of such func- 
tions determine the collective modes of a system, and 
within the TDHFA the effects of thermally-generated sin- 
gle particle excitations on these modes are included. The 
TDHFA has thus become a standard tool in condensed 
matter, nuclear, and atomic physics. 

Surprisingly, the TDHFA has not enjoyed such success 
in the area of Bose-condensed systems 2] . In part this is 
because there are other powerful methods for computing 
collective modes, based on the Bogoliubov approximation 
in which there is coherence among states with different 
particle numbers; i.e., where gauge symmetry is sponta- 
neously broken. Field theoretic techniques that incorpo- 
rate this were developed long ago, but in light of the 
significant progress on atomic gases in the last decade 
|^y, there is much renewed effort to improve upon these 
methods p|. Among these are generalizations of vari- 
ous mean-field schemes to properly include the effects of 
symmetry-breaking [||, and a "dielectric approach" that 
carefully treats Green's functions and response functions 
in a consistent way 0, ■ The challenge in these studies 
has been to treat the dynamics of the condensate and of 
the thermally excited (i.e., depleted) particles in a fully 
consistent manner. The resulting schemes are consid- 
erably more involved than the original Bogoliubov ap- 
proach. (For a review, see Ref.i) 

As we show below, one way to avoid the difficulties 
of disentangling the single-particle and collective excita- 
tions is to work with an ensemble in which the gauge 
symmetry is not broken |2j . In situations for which the 
particle number is relatively small, this may actually be 
preferable to broken symmetry approaches for which the 
associated fluctuations in the particle number can be- 
come problematic. Our chosen approach is the TDHFA. 



Its principle advantage is that the condensate wavefunc- 
tion and the single-particle excitations appear as states of 
a single time-dependent, non-local Hamiltonian, so that 
they may be treated on an equal footing with exchange 
effects fully included. The connection through a single 
Hamiltonian guarantees that conservation laws will be 
respected 0] , and we will demonstrate that the method 
correctly produces the gapless superfluid mode for a ho- 
mogenous system. 

Our principal conclusions are as follows. (I) We find 
that in order to properly deal with exchange, one must 
adopt a constrained grand-canonical ensemble for the 
density matrix first introduced by Huse and Siggia [Tlj . 
in which the number of particles in the condensate is fixed 
while the occupation probability for other single particle 
states is given by the standard grand canonical ensem- 
ble. Within this ensemble, one finds a gap in the single 
particle spectrum between the condensate and the other 
single particle states 0, > the latter of which we call 
depletion states 0|. This exchange gap is analogous to 
the single particle gap that arises in superconducting sys- 
tems and only arises when the system is Bose-condensed, 
and thus may be viewed as a superfluid order parameter. 
The gap should be observable in a tunneling experiment 
between a normal Bose gas and a Bose condensate. (2) 
We define a depletion Green's function G(q, to), whose 
poles occur at the single-particle energies of the deple- 
tion states in the Hartree-Fock approximation, which by 
construction does not have a pole at the collective mode 
frequency, as is expected for the full Green's function. 
The response of this Green's function to a weak, time- 
dependent perturbation, coupled with the equation of 
motion for the condensate wavefunction - which is a finite 
temperature generalization of the GP equation - allows 
us to generate equations of motion for the condensate and 
depletion states governed by the same effective Hamilto- 
nian, and to define response functions whose poles occur 
at the collective excitations of the system. This natu- 
rally captures the interplay between the condensate and 
depletion states. (3) We solve these coupled equations 
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for the simplest case of an infinite uniform Bose gas, and 
demonstrate that whenever there is Bose condensation, 
there is a gapless collective (zero sound) mode. This is 
usually identified as the superfluid mode in approaches 
where the gauge symmetry is broken; in our approach we 
find the mode even though the symmetry is kept intact. 
We demonstrate that the density response function con- 
tains structure that could not arise in simple Bogoliubov 
approaches where only the collective mode is retained. 

This paper is organized as followings: In Sec. [HI we 
derive the static and time-dependent HF equations at fi- 
nite temperature by using a constrained grand canonical 
ensemble. We then show how this can be used to define 
the depletion Green's function which captures the ener- 
getics of particles outside the condensate. In Sec. Mil 
we demonstrate how the depletion Green's function may 
be used to generate an approximate form for response 
functions. We also show that the same result may be 
obtained directly from wavefunctions, firmly establishing 
the connection with the TDFHA. In Sec. IIV1 we apply 
the general results in previous sections to a homogenous 
system with a contact two-body interaction. Discussions 
and conclusions are presented in the final section. Fur- 
ther details are presented in the Appendices. 



II. STATIC AND TIME-DEPENDENT 
HARTREE-FOCK EQUATIONS AT FINITE 
TEMPERATURE 

We begin with Bose particles in an external potential 
U(r) that is time-independent. The Hamiltonian is 



where 



H= drft(r)H ip(r) 



+ 2 j dr x dr 2 ^ {ri)^ {f 2 )V {fx - r 2 )i>(r 2 )'ip(r 1 ), (1) 

where Hq = — — h U(f) is the non-interacting Hamil- 
tonian and V{f\ — r 2 ) the two-body potential. For a neu- 
tral Bose gas, V is usually short-range and can be taken 
to have a contact form V(r) = gS(f) if the gas is dilute. 
We wish first to find eigenstates of H in the Hartree-Fock 
approximation (HFA) at finite temperature |lllll2|. 

We begin with the standard HFA, which we shall see is 
fine for high temperatures, but becomes a poor approx- 
imation when the system is Bose-condensed. We seek a 
single-particle Hamiltonian 



Hhf — / ,£af3a a a/3, 

a/3 



(2) 



where the indices a label states of a single-particle basis, 
which minimizes the free energy. In terms of these states 
the Hamiltonian Eq. JQ) can be written as 



H 



a/3 



alap 



E 

a/3 <57 



Vaps-^a^a^as^, (3) 



h 2 v 2 

2m 



U(r) 



M^) (4) 



is the matrix element of the non-interacting Hamiltonian, 
with ip a the single particle states, and 



V a/ 35-y = J dfi J df% ip* a (r*i )ipfi{r2)V{?i - t* 2 ) V><5 (^2)^7 (ri) 

(5) 

is the matrix element of the two-body potential. 

At finite temperature fcsT = 1/(3 where ks is the 
Boltzmann constant, the expectation value of an oper- 
ator O is 



(6) = TrDO, 



(6) 



where D is the exact density matrix operator. In the 
grand canonical ensemble this is 



D = 



(7) 



Tre -/3(H- M 7V) Z 

where Z is the corresponding partition function 

Z = Tre-K"-^l (8) 

The Hartree-Fock approximation is based on the varia- 
tional principle that for any trial density matrix D r , 
one always has for the free energy 



n[D r ] >n[D}. 



(9) 



In practice one chooses a form for D r that minimizes 
the free energy while allowing calculations with it to be 
tractable. In the Hartree-Fock approximation we choose 
the trial density matrix 



1 



,-0(H HF -^N) 



with 



(10) 



(11) 



The parameters e a p and _ib a are determined by minimiz- 
ing the trial free energy [l| 

^var — ~~q ln^^jar r ^- 1 ^-^) var H var -\- TrZ) i)ar _ff. (12) 

Defining the single-particle density matrix p a p as 

Paj3 = ^D var a) a a,3, (13) 

one has 

TvD var H HF = Tips. (14) 
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The variation of Q var yields 

6tt var = - — Tre - ^""- 6H var - TrSpe - TrpSs + 6(H) 

^var 

= -Tre6 p + 6(H) =Tr I ^p--e) 6 p, (15) 



with 



( H ) - ^2 £ al3P a + ^ E ^"/37<5 [PayPf36 + Pa6Pf3~/} ■ 



a/3 



a/3 (57 



(16) 



By requiring 6£l var — 0, one arrives at 
d(H) 



dp/3a 



= £ a/3 + VaiPSPSi, (17) 



7<5 



where V ai p$ = V ai p$ + V a7 sp- This is the Hartree-Fock 
equation in an arbitrary single-particle basis. 



A. Constrained Grand Canonical Ensemble 

The above derivation is valid for all temperatures, and 
the condensate plays no special role. However, the re- 
sult has significant problems for a Bose-condensed sys- 
tem: (1) It predicts macroscopically large fluctuations in 
the total particle number of the system below the critical 
temperature, which is unphysical. In the diagonal basis 
where n a = p aa = (e^ £< *~' 1 )) -1 , the source of this prob- 
lem may be traced to the parameter e Q= o — P associated 
with the condensate wavefunction, which becomes arbi- 
trarily small in the thermodynamic limit, (ii) With a con- 
tact interaction, the grand canonical ensemble produces 
the same mean-field potential for all particles, and does 
not yield the expected appearance of an exchange energy 
only when particles are in different states, not when they 
are in the same state. This is not a serious problem when 
the system is above the critical temperature since each 
level is microscopically occupied. However, this is a poor 
approximation when the system is Bose condensed. 

To correct these problems, one may introduce a con- 
strained trial grand canonical ensemble [Tl| 

D' = -^—Tre- ^^ s «-^< aa 6 t w . (18) 

var yi a^aojNo v ' 

var 

with Z' var chosen as usual to normalize the distribution. 
In this expression we have expressed Hjjp in a diagonal 
basis. This ensemble essentially excludes the dangerous 
condensed mode from statistical averaging , since this 
is what causes the problem in standard HF. As we will 
see, this removes the spurious exchange energy among 
particles occupying the condensate mode, while keeping 
it among particles occupying different levels. Consistent 
treatment of exchange in a Bose condensed system turns 



out to be essential for obtaining the expected gapless 
superfluid mode of an infinite uniform system. 

Using Eq. I|18|l . one gets the variational free energy 



ttvar = -pN + ^ [kBTln(n a + 1) - e a n a ] + (H), 



a#0 



(19) 



where (H) is 



1 x - 



i. 



(H) = ^s° aa n a + - ^ V a p a pn a np - -Vbooo-ZVo<5a,o 



and 



a/3 



n a = N 8 a n + 



(20) 

6 a -£Q. (21) 



When Q var is minimized with respect to Nq and n a (a ^= 
0), which is equivalent to variation with respect to e a , one 
finds 



(22) 



dln(n a + 1) de a d(H) 

° = ^ e a -n a — + —, a^O. (23) 

These equations are straightforwardly solved by setting 
e a = e° aa + V a p a p(l - 6pfl6 a p)np + Voooo-Nb<S 0) o 

= e° aa + S QQ , for all a (24) 



and 



P = £o, 



(25) 



where 



£ QQ — ^""^[Vaffaff + V a pf} a ]np — Voooo-Wo^a,o (26) 



is the diagonal matrix element of the Hartree-Fock self- 
energy. 

The problem is not fully solved because we do not have 
expressions for the wavefunctions that determine the ma- 
trix elements V a p^s- To find these, we minimize the free 
energy with respect to the single particle wavefunctions, 
keeping in mind that they must form a complete or- 
thonormal set |ll| . This constraint may be enforced if 
we write the variation in the form 

6lp a = ^ Va-flp-y, Or 6a^ a = ^ Varf a \i ( 27 ) 



where = —77*^. The resulting variation of the free 
energy £l var may be written as 



6Q, var = r) VfJ ,([H,a'la u ]). 



(28) 
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The proof of this result is given in Appendix [XJ and is Eq. (|31|l suggests that the self-energy matrix should be 
valid whether or not the system is Bose condensed. Now defined as 
substituting Eq. J3J) into Eq. (J2EJ, one gets 
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(29) 



The free energy should be stationary with respect to vari- 
ations in rj, 



= o, 



which results in the equation 



(30) 



(32) 



which we note is consistent with our earlier definition of 
the diagonal self-energy matrix elements. In real space 

(31) 

this may be written as 



£(fi,f 2 ) = ^2 (Tl)^^) = 5(n - r 2 ) dr 3 V{r 2 - r 3 )p(r 3 ,r 3 ) + V(fi - r 2 )p(ri,r 2 ) 



N ldr 3 W (r 3 ,r 3 ) [V{r x - r 3 ) + V(r 2 - f 3 )] W (n,r 2 ) + N IfoooWbft ,r 2 ) (33) 

I 



where p(f\, r 2 ) is 

P{ri,f 2 ) = ^2n a ip a (r 1 )ip^(f 2 ), (34) 

a 

and W a (ri,f 2 ) is defined as 

W a (r 1 ,r 2 )=i> a (r 1 )i>* a (r 2 ). (35) 

With these definitions, and using the completeness of the 
basis, Eq. (|31|l is equivalent to the Hartree-Fock equation 



V 2 



i> a (r) + ld7Z(?,r*W a (7) = e a ip a (r) 

(36) 



with 



£ a = £aa + X] ~ -^O^OOO^a.O- (37) 



Note this is consistent with Eq. (|17|l expressed in a diag- 
onal basis. 

The form of the self-energy in Eq. (|33|l . as most clearly 
expressed in Eq. I|31|l . has an interesting consequence: 



the self-energy for the bosons in the condensate is dif- 
ferent than that of higher energy states. For the case 
of a uniform, infinite system this means there is a gap 
between the single-particle energy of the condensate and 
those of the excited states. This property of Bose con- 
densates is well-known |l5| . and finding an appropriate 
way to deal with it is one of the major challenges in de- 
veloping approximations for the excitation spectrum of a 
BEC 51. For the case of a uniform infinite system, this 
gap is present in the single-particle spectrum in spite of 
the expected gapless collective mode spectrum. In order 
to deal with this, it is helpful to develop different Green's 
functions which capture one or the other portion of the 
excitation spectrum, as we now proceed to do. 



B. Depletion Green's Function and the TDHFA 

In formulating a TDHFA, it is useful to define Green's 
functions in imaginary time and consider self-consistent 
approximations to their equations of motion |l0j . In the 
fermion case, poles of the Green's function in the ab- 
sence of a time dependent potential give the spectrum of 
number-changing excitations, while the response of these 
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Green's functions to time-dependent potentials give col- 
lective excitations. This allows one to conveniently sep- 
arate out these sectors of the energy spectrum. In a 
Bose condensed state, these sectors become entangled 
in the standard Green's function because one may add 
a particle to the condensate and then excite a collective 
mode, yielding poles at collective mode frequencies. Dis- 
entangling the single-particle spectrum from the collec- 
tive mode spectrum in the Green's function then becomes 
quite challenging. 

As we now demonstrate, the TDHFA for Bose con- 
densates can be developed in a way that is analogous to 
what works so well for fermion systems. To to this we 



define a depletion Green's function, incorporating all the 
information about the single particle states other than 
that of the condensate. Within the static Hartree-Fock 
approximation this has the form 



G(fi,f 2 ;iw„) 



iu! n - £ a + P 

a^tO 



(38) 



Writing this Green's function in imaginary time, 
G(y 1 t 1 ,y 2 t 2 ) = ^X)n e ~^" (ri_T2)( 5'(ri,r2;iw„) one may 
easily show that it satisfies the equation of motion 



ffo(fi) 



Gfar\,f 2 T 2 ) - dr T 1 (ri,r)G(rT 1 ,f 2 T2) = 8 fa - T 2 )[Sfa - r 2 ) - ipo(fi)ipo(r 2 )]. 



(39) 



By excluding the condensate state from the depletion Green's function, we avoid the process that leads to poles at 
the collective mode frequencies. The lowest energy poles then reflect the single-particle spectrum. 

To compute collective modes of the system, it is convenient to look at Green's functions in the presence of a 
time-dependent potential [T^j . The natural generalization of Eq. (|39|l to this situation is 

— -Hoif^-SUifun) 



G(fiTi,r 2 T 2 ) - jdr /drE(fm, rr)G(rT, f 2 r 2 ) = 8fa - T 2 )[8(fi - r 2 ) - tpo(fi,Ti)ipQ(f 2 ,T 2 )}, 

(40) 



where SU is a time-dependent potential which we will ultimately treat perturbatively. Note that in writing down this 
equation, the self-energy Eq. Q33J) now has time dependence, and is explicitly given by 



Z(riT 1 ,r 2 T 2 ;8U) = 5 fa - t 2 ) 



8(n ~ r 2 ) dr 3 V(r 2 - r 3 )p(r 3 T ll r 3 T l \8U) + Vfa - r 2 )p(r*iTx,r 2 Ti;8U) 



8 fa - t 2 )N \ dr 3 [Vfa - r 3 ) + V(f 2 - r 3 )] Wo^n, f 3 n; 8U) - V ooo > W (fin,f 2 n;8U) 



(41) 

The time dependence enters through the wavef unctions in the quantities p and W a (Eqs. i|34|l and (|35|l 'l. and we have 
noted that these quantities are now functional of 8U which is ultimately responsible for the time dependence. 

It is useful to note at this point that we have made a crucial assumption, which can be understood as the essential 
underlying approximation of the TDHFA: we allow only the wavefunctions to change with time, while the occupations 
n a remain stationary and equal to their values for SU = 0. This can be shown Q to be equivalent to an assumption 
that the entropy of the system remains unchanged in the presence of SU . 

Noting that Eq. 1341 may be recast in the form 



p(fiTi,f 2 T 2 ) = G(r 2 ,r 1 + ,ri,Ti) + N W (fiTi,f 2 Ti) 8fa-T 2 ), 



(42) 



with Wo(riTi, r 2 T\) = ipo(riTi)ipQ (r 2 r 2 ), we see that Eqs. (|33|) . (|40(l and 14211 nearly form a closed set of equations. 
We have left to determine the time dependence of tpo- An important aspect of the problem is to assure that our 
TDHFA obeys particle conservation, a feature that is often difficult to build into collective mode calculations for Bose 
condensates In the present case we can guarantee this by having the condensate wavefunction controlled by the 
same effective Hamiltonian as the excited states, via Eq. (|40() . Thus we take 



~-H {?i)-SU(fi,n) 



it'ofan) - Jdr JdrJ:(fiTi,r,T)ipo(fT) = 0. 



(43) 



With this equation for Vo, it is easy to ver- ify that the overlaps J dr^Q(rr)G(rr, r'r') and 
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J dr'G(rT,r'T')ipo(r'T') vanish, so that the deple- 
tion Green's function involves no change in the number 
of condensate particles even in the presence of the time- 
dependent potential. It is the possibility of changing 
this and simultaneously creating a collective excitation 
that allows the collective mode spectrum to appear in 
the standard Green's function. Thus by working with 
the depletion Green's function we avoid the entangling 
of particle-conserving excitations and single-particle 
excited states that characterize the approaches based on 



broken gauge symmetry. 

In principle these equations may be solved self- 
consistently to develop a mean-field approximation for 
this many-body system in a time-dependent potential. In 
practice this is a formidable task, so one instead focuses 
on the linear response of G and ipo to small perturbations 
SU. These may be used to construct, for example, the 
density response function, whose poles give the collective 
modes of the system [Toj . 



I 

III. LINEAR RESPONSE 



A. General Formulation 



We begin by expanding Eq. 1401 for small SU, retaining only terms that are linear in this quantity. This leads to the 
equation 



SU '(fiTi)W 'o(ri,r 2 )8{T\ - r 2 ) 
d 



H H f{t2) 

OT 2 



dn 



dr 2 



HHF(r 2 ) 



S (n - t 2 )S W (n , r 2 ; n ) 



SG(fi,T 1 ;r 2 ,T2) - df f a S^(fi,ri;f 2 ,T 2 )Wo(f 3 ,f 2 )5(Ti - t 2 ), 



(44) 



where Wo(ri,r 2 ) — 8(ri — f 2 ) — Wo(fx,f 2 ), and SWo(ri,r 2 ;r) = Stpo(ri, t)V>o (^2) + ^0(^1)^0(^2,1"). The operator 
Hhf(t) has the meaning, for example, 

HHF{r2)5Wo{r u r 2 ;Tx) = H {r 2 )6W {r u r 2 ; n) + fdrS(f 2 ,r)SW (n,r;n) 

I 



with £ given by Eq. (|33l) . The variation of the self-energy, 
<5E, comes from the fact that the wavefunctions ipa are 
functional of SU; we will provide an explicit expression 
for the specific case of a contact potential below. Be- 
fore proceeding with this, we demonstrate that the ideas 
developed above may be used to compute an important 
quantity, the density response function. 

B. Density Response Function 

An alternative procedure for avoiding the singularity in 
the standard imaginary time Green's function is to work 



directly with the wavefunctions and the density matrix 



a 



which is a generalization of Eq. (|34l) . In this context it 
is convenient to work with real rather than imaginary 
time. The equation of motion for real time density ma- 
trix p(r*i, r*2, t; SU) can be easily obtained from the time- 
dependent Hartree-Fock equations for the wavefunctions, 



at 1 



-p- + U(f 1 )+SU(f 1 t 1 ) 



with the result 



ipa{fiti) + dr 2 dt 2 Z(riti,r 2 t 2 ;SU)'ilj a (r 2 t 2 ), 



p(ri,r 2 ,t;SU) 

dr 3 [Efti, r 3 t; SU)p(f 3 , f 2 ,t; SU) - p(r u r 3 , t; SU)Z(r 3 t, r 2 t; SU)} 



(45) 



" d 








\2m 


2m J 



(46) 
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In writing this equation we have analytically continued 
the self-energy to real time, which in practice simply in- 
volves replacing r — ► t in all the arguments of Eq. (|41|l . 

Expanding p{r\, r 2 , t; 8U), T,(fit,r2t;5U) and 
Wo(fit,f 2 t;8U) around their static HF values and 



retaining only terms that are first order in SU(ft), one 
obtains a linearized equation for Sp, 



[SUif^-dUi^jpin,^) = 



d - 
ig^ - H HF (ri) + H HF (f 2 ) 



5p(nt,r 2 t + ;U) 

'dr 3 [Spirit, r 3 t;SU)p{r 3 ,r 2 ) - p(fi, r 3 )SJ:(r 3 t, r 2 t; 8U)] 
A density matrix response function 10] may be defined as 



X H (nr 3 ,r 2 f 3 ,t~t 3 ) =6(t-t 3 ) 



Sp(fi,f2,t; SU) 



SU(r 3 t 3 



(47) 



(48) 



su=o 



with the more standard density-density response function then given by X^fiTs, fxf%, t — t 3 ). One can rewrite Eq. 
(B7|l as 



d - 
ig^ - H HF (ri) + H HF (r 2 ) 



X R (r 4 r 4 ,r 2 r 4 ,t- t 4 ) = [S(fi - r 4 ) - S(r 2 - f 4 )} 5(t - t 4 )p(ri,r 2 ) 



where 



r K (fif 4 , r 2 f 4 , t-U) = 9(t - U) 



+ jdr 3 [r fl (rir 4 ,r 3 r 4 ,i-t 4 )p(r3,r2)-p(fi,f3)r fl (r 3 r 4 ,f2,r ; 4 ,t-i 4 )] , (49) 
5T,{r l t 1 r 2 t 1 8U) 



8U{r 4 t 4 



su=o 



S(n - r 2 ) fdr 3 V(r 2 - f 3 )x R (f 3 f 4 ,f 3 r 4 ,t - t 4 ) + V(fi - r 2 )x R (rir 4 , r 2 r 4 , t - t 4 ) 

SW (r 3 t,f 3 t;SU) 



N df 3 [V(fi - f 3 ) + V{f 2 - f 3 ) 



5U(f 4 t 4 ) 

N>{ /'//". Ttr, • /".. ! ~-V{r 2 -r 3 )}W (r 3 ,r 3 )-Voooo 

6W (r 3 t,r 3 t;SU) 



W (n,r 2 ) 



8U=0 



+ 2N df 3 dr 5 V{r 3 -f 5 ) 



5U(f 4 t 4 ) 



8U{r 4 t 4 ) 
Wo(f 5 ,r 5 )W Q (ri,r 2 ). 



su=o 



(50) 



su=o 



The functional derivative of Wq with respect to 8U{rt) requires one to work directly with the wavefunctions. This is 
discussed in detail in Appendix B, with the result (Eq. (|B10(l 'l 



5W Q {r 1 t 1 ,f' l t 1 ;8U) 



SU(r 2 t 2 ) 



5(7=0 



-E 
+ E 



2n lu + Ep - so 
2ir uj + Eq — ep 



[M/»a + r£ 0)A/ »] r^M^MnWoi^). (5i) 



To make further progress, it is convenient to expand the functions x R an d F R in terms of the static HF eigenstates 

ip a as 



/oc 
2^Xe*/3,\tx.{u) 



a/3 
A// 



/oc T 



iw(tt— i 3 ) 



(52) 



a/3 
A// 
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Substituting these expansions into Eqs. (|49|l . (|50|l . and (|51H . one obtains 

+ e M - ex) Xxwpfa) = ( n v ~ ™a) + rf MiQ/3 (w)] 

and 



(53) 



Ai^i 



Ai/ii 



A1M1 



Ai^ii 

+ A/bVoooo ($\o - ^o) Xa^/jH- 



(54) 



Using Eq. (|54|l in Eq. (|53[) to eliminate T, one ultimately may write what amounts to a matrix equation for X\n «a(' 4 ')) 
(w + e M - £a)xv :Q/3 (w) = (rty - nx)S a x5p ll + [n„ ~ n\) Y [Vx^Xm + ^a M imaJ Xa iW ,^( w ) 



Ai/ii 



(n p - 7l A ) [VViA ljU (<*Ai,0<*A,Ai + <W<W) + ^A Mlj uAi (<*A,0 + ^,o) (^A^A + <W/J] Xa iMi ,a/3 (<*>)• 

Ai/ji 

(55) 



r 



Interestingly, the first line this equation has the same 
form as what is obtained from the TDHFA using the 
standard grand canonical ensemble, although the ener- 
gies and occupations are different because of the non- 
local terms in the self-energy that arise from using the 
constrained grand-canonical ensemble. The terms in the 
second line also come from these non-local terms. This 
equation for the density response function provides a nat- 
ural way to describe the coupling between the conden- 
sate and normal components through the matrix indices, 
which may refer to the condensate state or to the de- 
pletion states. We will illustrate this explicitly when we 
solve Eq. (|55|) for a homogenous system below. 

Finally, it is useful to note that Eq. (|55|) can written 
in the form of a Bethe-Salpeter equation 

XAm^M = XA/^a/jM 

+ XA^pr/M^^HX^a/sM, (56) 



where a| g(w) is defined as 

X%, a p(u) = 
and the kernel K as 



U + IT] ■ 



£A 



SaxSpfi: (57) 



+ Vpau n [1 - (Spa + M (8 V p + 5 au )] (58) 
which is independent of the frequency. 



Equation (|55|) is valid for a dilute Bose gas with any 
shape of static external potential U (r) and a general two- 
body interaction potential V(f\ —?2j. In general it must 
be solved numerically. In the next section, we will show 
that it can be solved analytically for a uniform, homo- 
geneous system of bosons interacting via a contact two- 
body interaction. 



IV. APPLICATION TO A HOMOGENOUS 
SYSTEM 



In this section, we apply the above results to dis- 
cuss the single-particle and collective excitations of a ho- 
mogeneous system with a contact two-body interaction, 
V(f\ — r%) — g8(fx — r*2). This is a common and quanti- 
tatively accurate approximation for dilute Bose gases in 
many situations. We begin by reviewing the solutions of 
the static Hartree-Fock equations for this case Q . 



A. Static Hartree-Fock for a Bose-Einstein 
Condensate 

For a uniform homogeneous system (U(r) — 0), 
the Hartree-Fock Hamiltonian is diagonalized by plane 

waves, ip%{r) = ~^ elk ' r > with ^ the volume of the sys- 
tem, and the single-particle energies are easily shown to 
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be 



where 



2 9P - gpoh o 



(59) 



where = p = N/Cl and p = N /fl. Below the 

critical temperature, the chemical potential is given by 
the lowest eigenvalue 



p = s a = g(2p - p ). 



(60) 



Note that there is a gap in the single-particle spectrum 
A = gpo with respect to the chemical potential for a Bose 
condensed system. As discussed in the Introduction, if 
the grand canonical ensemble is used, there is no such 
gap and the corresponding single-particle energy is = 

+ 2gp. However the absence of the gap is an artifact of 
strong number fluctuations in the condensed state, which 
are more properly controlled by the constrained grand 
canonical ensemble. 

The gapped structure of the single-particle spectrum 
may be probed in principle via a tunneling experiment. 
One interesting possible approach involves a "single-atom 
pipette" that has recently been proposed 0, in which 
bosons are loaded into a strongly localized potential. Be- 
cause of the potential well, there can be a considerable 
gap for excitations within the pipette, and at low tem- 
peratures all the atoms will be at the same energy. By 
bringing the pipette sufficiently close to the bulk BEC for 
a fixed time, atoms may tunnel out of the pipette into the 
bulk BEC, and a measurement of the number of atoms 
remaining in the pipette after this process allows one to 
infer the tunneling rate. For a relatively deep pipette 
potential, one expects the atoms to tunnel only into the 
condensate state of the BEC. As the energy is raised, 
however, there should be a sharp threshold energy, set 
by the gap, above which atoms may also tunnel into the 
excited states. This will show up as a non-analytic contri- 
bution to the tunneling rate reflecting the density of de- 
pletion states in the BEC. For a bulk, three-dimensional 
BEC, this would be proportional to (E — Eth) 1 ^ 2 , with 
E the energy of the bosons in the pipette, and E t h set 
by the gap energy. We note that the exponent might be 
renormalized by shakeup effects |17| . but we expect the 
non-analytic behavior to be robust and to reflect the gap 
in the single-article spectrum. 

The critical temperature temperature T c is determined 
by taking the limit p Q — ► 0, so that the number of par- 
ticles depleted from the condensate tends to the total 
number of particles N. This leads to the condition 



N 



V 



which gives 



k B T c 



(2tt) 3 



47T 

2m' 



d A k- 



1 



Mi 



,2/3 



33/2(1) 



2/3 



07 ( X ) = Yl 



71=1 



(63) 
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FIG. 1: (Color online) Condensation fraction as a function of 
the temperature at two dilute parameters. 

The condensate density po below T c is found by re- 
quiring depletion of the condensate to be equal to the 
thermal occupation of the excited states Q, 



1-* = 



dk- 



= T73ff3/2(^), ( 64 ) 



p (27rfpJ a e Ml+9Po)_ 1 \ T 
where z = e~^ 9Po , and 

X- =(__f) V2 „./3 = W / a (65 ) 
\ m J 

is the product of the thermal wavelength and the inverse 
average distance between particles. We note that al- 
though the Hartree-Fock approximation yields precisely 
the same T c as for a non-interacting system, the phase 
transition into the condensed state is first-order rather 
than continuous. This can be seen in the behavior of the 
condensation fraction as a function of the temperature, 
as illustrated in Fig. ^ 

From Eq. I|64|) and Fig.^ we can see that in the HFA, 
the repulsive interaction enhances the condensation rel- 
ative to the non- interacting case [l5| . 



(61) B. Collective Modes from the Depletion Green's 

Function 

We would like to solve Eq. (|44|) for the case at hand; 
in order to do this we need to specify 5E. For contact 

(62) interactions in the homogeneous system, this is easily 
evaluated, with the result 
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<5£(fi, ti; r 2 , t 2 ) = I 2g [8p(fi,Ti) + 5po(fi,Ti)] S(n - r 2 ) + g 



5w(fi,f 2 ;Ti)p (fi,f 2 ) - -^Sp (fi,f 2 ;Ti) > 5(n ~ t 2 ). 

(66) 

In this equation, Sp(fi,Ti) — 5G(riri, fir* ) is the variation of the depleted particle density, Spo(fi,Ti) = 
NoSWo(fi,ri,Ti) is the variation of the condensate density, and Sw(fi,f 2 ;ri) = —[Sipo(ri) + <5-0o(^i) + if>o(r 2 ) + 
5ipQ{r 2 )]/ As in the case of the density response function, we can make progress by expanding the various 
quantities in terms of the unperturbed single particle states. Thus we dehne 8G^^{iuJi,iL0 2 ) via 



6G(r lTl ,r 2 T 2 ) = - £ e --^+^_ £ ^-^^^(i^ia*). 



(67) 



It is useful to notice that the Fourier transform of the depletion density, 6p(k,iu) n ) — Jq dr J dfe l ^' r+zlJJnT S p(r, t) , 
satisfies 



Sp(kiw n ) = — SG k+k'M'( iuJn + iuj 'n> iuj n)- 



(68) 



Using Eq. (|67(l . Eq. (|44|l may be recast in a particularly simple form. For k\,k 2 ^ 0, we find 
[e j j i -twi][e fe -twa]<yGjj ii& (iwi > tW2) + ^ ^ *Gg, +Ajfij j>(tw / +iAa;,<w / )+2ff*/Jo(Afc > iAw) = -5U(Ak,iAu), (69) 

where Acu = cui — u)2, Ak = k\ — &2, and 



6p {Ak,iAcj n )=-j dr J dre lAk -r +lA ^ T {6Mr,' 
SU(Ak,iAu; n ) =i / dT / dre^-^+^-^C/^T; 



fi./o 



It is clear that in Eq. (|69|l we would like to write k\ = k 2 + Ak, u)\ = Alu + uj 2 , divide by the equation by [g£ — 
iwi][e^ — itu 2 ] and sum with respect to k\ and wx- However, there is a caveat: in doing this our indices run over 
&x=0 and k 2 = 0, which introduces further terms in Eq. (|69|l . It is not difficult to show that these terms vanish in 
the thermodynamic limit (fl — > oo), which is the case of interest to us. Performing the above steps, Eq. (|69[) becomes 
a simple linear equation, 



P(Ak,iAuj) 



25 



5p(Ak,iAw) - 2gSp (Ak 1 iAuj) = SU(Ak,iAoj), 



(70) 



where P(q,iu n ) = ^ _n s(%)]/[iw„+£ : jj + --eg], eg = eg-^t, n s (s) = l/(e /3e -l)is the Bose occupation 

factor, and the prime on the sum indicates that k, k + q = should not be included. 

Since the problem involves two density disturbances, the depletion density Sp and the condensate density Spo, we 
need a second equation. This can be obtained directly from the groundstate wavefunction disturbance, as described 
in Appendix B, since Spo(r, t) — [5^o(f f , t) + St/j^f, r)]. Using the method of Appendix B, one may easily show 
for imaginary time, 



1 



ft 3 / 2 



5U(k, iuj n ) + ST,(k, iu> n ) 



lUJ n 



(71) 



where the prime on the sum indicates k = should not be included, and 

5E{k,iu n )= [ dr x dr 2 [ dT 1 dT 2 e iu " T1+i ^S'E(f 1 T 1 ,f 2 T 2 ), 



(72) 
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Substituting Eq. ® in to Eq. J7TJ, 

one obtains 
1 



VI! V Si 



where Stp (k,iuj n ) = J dff dre lk ' T '5^p (r,T) and StpQ(k,iuj n ) = J dfj dTe lk " r Stp^r, t). Similarly, one finds 



1 2g - 

—7=SU (k, ioj n ) + —=5p(k, ioj n ) + gp 5i/j Q (k, iui n ) 



Equations f?3*)) and (f74*j) may substituted into the definition of 5po to obtain 



£ k ~ 9PQ 



£ k - 9Po 



5p (k, iuj n ) - 4gp 5p(k, iu n ) = 2gp SU(k, iu n ), 



(73) 



(74) 



(75) 



providing the second equation needed to compute the density disturbance. 
Equations (|70|l and l|75|l may be combined into a single matrix equation, 



25 



-2<? 



- 2 .9 Wo 



etf—gpo 



Sp(q, iuj n ) 
Sp (q,iuj n ) 



5U(q,iu n ). 



(76) 



Collective modes of the system propagate when Sp, Spo ^ even if 5U = 0. This is only possible if the matrix on the 
left hand side of Eq. I|76(l has vanishing determinant, 



25 



Sq ~ 9Pa 



(Sq + gpo) 



- 8g 2 Po = 0. 



(77) 



Upon analytic continuation (iui n — > u> + iff) , Eq. i|77|) supports a linearly dispersing gapless mode (zero sound) for any 
pa =/= 0, as expected for an infinite uniform superfluid. This is a non-trivial check that our formalism obtains physically 
sensible results. Before analyzing this in further detail, wc demonstrate that the same result may be obtained directly 
from the linear response formalism. While somewhat more complex to carry through, this latter approach allows us 
to look at very general response functions, and so yields more information than the density responses above. 



C. Density-Density Response Function and Collective Excitations 



Our starting point for this analysis is Eq. (|55|l . For a homogenous system, \ R vanishes unless the indices are 
such that momentum conservation is respected. Thus we may set A = k\ — q/2, fi = k\ + <f/2, a = &2 — <f/2, and 
(3 = k-2 + q/2. Denoting 



(78) 



one obtains 



k\k 2 



gpo 



gpo (s iit _ f - ^ x£ l£a (#w). 



(79) 



Bringing the last term to the left side of the equation has particle spectrum, since 
the interesting effect of canceling the gap in the single- 



£ k 1+ i - e k t -l + 9Po ( 6 k u -i - 5 kt,i) - 4 1+ f _ ^-f 

(80) 
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where eS = is the free single-particle energy. This the self-energy used in the static Hartree-Fock analysis 
cancellation is only possible because our linear response E3- 
equation was generated in a way that is consistent with 



J 



Using Eq. J73}, one can obtain 



c 2 



; 0)2 



Nr, - n, r 



Nn - Tiff 



U k 2 ,y 



c 2 - ( £ £) 2 + 2gpoe% n ^ 



(81) 



Substituting this result back into Eq. (|79|) leads to 



= -7-70 TTo " 



()\2 



(4) 







7V - n g - 



fci+^g fci 



Now summing over k\ and fe, one gets 



c 2 - ( ff 0)2 + 2gP0 e° 



2ff 



P(?;w)- 



2fo4 



c 2 - (el) 2 + 2gpoe% 



(83) 



One finally may express x^„(<f;w) in a symmetric form 

P(q;uj) + P c (q;u) 



l-2g 
where we have denoted 



P(q;uj) + P c {q;u) 



(84) 



^ - (e") 2 + 2 gPo e° 



: 
9 



(85) 



which can be interpreted as the polarization function for 
the condensate. The poles of Xnn determine the collective 
modes, and are given by solutions to 



1 = 2g 



P{q;uj) + P c (q-iu) 



(86) 



It is not difficult to show that this is identical to Eq. (|77|l . 



D. Discussion 

Several comments about the linear response results as 
obtained from our TDHFA approach are in order. Firstly, 
At T = 0, P(q;uj) = 0, and the density-density response 



function becomes 



Xnn(q;uj,T = 0) 



with the pole given by 



l-2gPM"Y 



oj"(q ) T = 0) = 2gpe}+(s ,) 2 . 



(87) 



This is exactly the same as that obtained by linearizing 
the time-dependent GP equation. 

Secondly, we note that Eq. I|8t)|) yields a propagat- 
ing gapless mode for any po ^ 0. This means our ap- 
proach correctly captures the superfluid mode of the sys- 
tem whenever it is Bose-condensed. Thus in retaining 
the correct (gapped) structure for the single-particle ex- 
citations, we see that the gapless superfluid mode is not 
sacrificed. At T > T c , P c {q\uj) = 0, and writing P = P, 
the density-density response function is 



Xmi(g;w) 



with poles determined by 



P(q;u) 



l-2 5 P(<f;u;)' 



l = 2 g P(q;u), 



(89) 



(90) 
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FIG. 2: (Color online)Imaginary part of the density response 
for different values of q as a function of u> for a Bose conden- 
sate at inverse temperature j3 (in units of the inverse crit- 
ical temperature (3 C ) with density p and scattering length 
a = m/ATvh 2 g. 



which is the known result for the normal Bose gas |l8j . 

Thirdly, because of the structure of Eqs. lf76*|) and (|SrJ|l . 
one might expect to find two collective modes, one in 
which the condensate and depleted particle densities os- 
cillate in phase, and the other in which they are out 
of phase. These would correspond to zero and second 
sound, respectively. We find however that for attainable 
values of fi(T) and po(T) within the HFA, the two so- 
lutions of Eq. l|86|) occur with one at positive cu 2 and 
the other at a negative value; the second sound mode 
is thus overdamped. This is similar to the weak cou- 
pling limit [2|. However, unlike the weak coupling limit, 
here the second sound is overdamped because of the gap 
in the single-particle spectrum, which makes the deple- 
tion particle density too small to support a propagating 
second sound mode. The gap has another interesting 
effect on the interaction between the condensate and de- 
pleted particles. This is illustrated in Fig. |2 which de- 
picts the imaginary part of the density response function, 
Xnra(<7) w + it]) f° r fixed q as a function of u>. Two fea- 
tures are prominently visible: a peak (which becomes a 
delta function in the limit q — > 0) representing the super- 
fluid mode, and a zero in the response that appears for 
q 2 /2m > 2A = 2gpo. The width of the superfluid mode 
arises because of the interaction of the condensate with 
the depletion particles, and vanishes at low temperature 
as e - A / fc sT Note we can obtain this temperature depen- 
dence only because our choice of self-energy creates a gap 
between the condensate and excited states, even though 
we use a single Hamiltonian to describe them. The zero 
arises when J 1 — (^j) (^- — 2 A), and represents a fre- 
quency at which 5p in Eci. 1761 vanishes. In this situation 
the condensate precisely screens the perturbing potential 
8U for the depletion particles. Since the depleted parti- 
cles are unperturbed and the condensate cannot absorb 
energy away from the superfluid frequency, no energy can 
be absorbed by the system, leading to the zero. It is in- 
teresting to note that an observation of this effect would 
allow one to measure the energy gap of the system. 

Finally, for small depletions p, the velocity of the zero 
sound mode found from Eq. I|86|l may be shifted either 



upward or downward, depending on the numerical value 
of the gas parameter. This is most easily demonstrated 
by expanding the equation for small p and small P. The 
correction to the sound velocity can then be shown to 
have the form Sc = j^{gAp + 8g 2 pP }, where c 
is the zero sound velocity at zero temperature, Apo the 
change in the condensate from its zero temperature value, 
and Pq = lim 9 _>o P q (u> = coq). For small values of p, 
Pq < 0, and the mode velocity decreases with tempera- 
ture, as is commonly found in the Bogoliubov approxima- 
tion. By contrast, we find for larger values Pq > 0, and 
for large enough p, its value is sufficiently large to render 
Scq > 0. This means relatively dense superfluids may 
have increasing zero sound velocity with temperature, as 
is found in the RPA and appears to be consistent with 
data for 4 He [l9|. That the TDHFA can capture both 
these limiting behaviors demonstrates the utility of the 
method. 



V. CONCLUSION 

In this work, we demonstrated that the TDHFA may 
be developed in a way that does not break the gauge 
symmetry of the Hamiltonian, thus allowing us to cor- 
rectly obtain a gapped single particle spectrum, and yet 
correctly produces a gapless superfluid mode in the col- 
lective response of an infinite uniform system. By devel- 
oping both the depletion Green's function and response 
functions, we can examine the single-particle and collec- 
tive mode spectra separately. The key to obtaining these 
features in a consistent way was to retain the nonlocal 
terms in the self-energy (see Eqs. (|33[) and J2J), which 
were required by the orthogonality of the single-particle 
wavef unctions, and ultimately lead to a precise cancel- 
lation of the gap in the single-particle energies when we 
calculate the density response function (Eq. iJSDJ). 

Finally, we point out that our equations for the infinite 
uniform system (Eq. I|76l) ) are formally similar to equa- 
tions for the density response obtained using ensembles 
in which the gauge symmetry is broken 0. This formal 
similarity however does not generally occur; for example 
if we had studied a uniform but finite size system (parti- 
cles in a box), we would find further terms that vanish in 
the thermodynamic limit, which do not appear in other 
approaches. Generally speaking, in computing collective 
modes of a Bose condensate, the challenge is to find how 
the condensate density couples to the depleted particles 
in a way that is conserving |10| , and preserves the gapless 
mode expected for an infinite system |5(. As discussed 
above, the TDHFA is guaranteed to be conserving as it 
is controlled by a single Hamiltonian (in contrast to many 
other methods), and we have demonstrated that it cor- 
rectly produces the expected superfluid mode. Beyond 
the case we studied in detail in presenting our method 
here, the TDHFA can be used to study inhomogeneous 
and/or finite size systems, and may be very naturally 
generalized to handle multicomponent Bose systems, ro- 
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tating systems, and even boson-fermion mixtures. Future 
studies will focus on these applications. 
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Now moving to the left side and noticing 

J"' | (at i ) na *~ 1 \0) = ^/n^~\n a , - 1) = a ai \n ai ), one ob- 
tains Eq. HA4f). 

The change of the free energy due to the basis change 

is 



APPENDIX A: PROOF OF EQ. lt28l) 

The off diagonal matrix elements of S can be found by 
varying the orthogonal basis tp a . We assume that I*) is 
a Hartree-Fock permanent constructed from orthonormal 
single-particle states ip a , which can be expressed as 



(Al) 



where |0) is a vacuum state. The variation of the single- 
particle states may be written as 

Sip a = ^77q 7 V' 7 , or 5a) a = ^?7 Q7 a 7 , (A2) 



where Ei is the unperturbed Hartree-Fock energy cor- 
responding to the Now substituting Eq. I|A4(I into 
above, one obtains 



Y,ie- 0{Ei -» N) {*i\[H,ala v ]\*i 



which is the result of Eq. . 



(A7) 



APPENDIX B: PERTURBATION THEORY FOR 
WAVEFUNCTIONS 



with the orthonormality of the basis ip a implying that 



El 



*7ot/3 + V*0a = 0. 



(A3) 



We can prove that the change of the permanent (|A1|I due 
to the basis change l|A2j) can be generally expressed as 



a/3 



(A4) 



This can be shown by substituting Eq. (|A2|) into Eq. 
(lAH . 

w=Ew n ^&«)--- i ^=^io>. 



(A5) 



In this Appendix, we seek a formal perturbation solu- 
tion of the TDHF equation since we need a relation of 
the functional derivative of W a with respect to SU. 

If the time-dependent external field 8U(rt) is weak, we 
can solve the time-dependent single-particle wavefunc- 
tions <p a (rt) using time-dependent perturbation theory. 
In order to do this, we expand the time-dependent single- 
particle wavefunction <j) a (ft) in terms of the static HF 
basis ip a (r) for times t > to — ► — oo as 



b a (rt;6U) = J2 C «p{t)e~ ie0t Mr)- 





(Bl) 



We seek a solution for C a p(t) up to the first order of 
5U(t). Substituting this into Eq. (|45|l . after expanding 
the self-energy E to first order of SU(ft) and using Eq. 
(I5H . one has 



Y,C a p{t!)e- ie e tl dr 2 dt 2 dr 3 

13 J 



6X^,^2, SU) 



J h SU(r 3 t 3 ) 



6U(r 3 t 3 )M?2). (B2) 



8U=0 
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We can rewrite Eq. (|B2|) as 

a J J — oo 



Multiplying both sides of Eq. (|B3jl by ijj* (ri ) , integrating over r\ , and using the Fourier expansion 

doj 



5U(ft)= l^-5U{r-uj)e-^\ 



one gets 



: d C a p{ti) _ v fdw 



e~, — e/3)ti 



Writing 



C a p(t) = S a p + C^(t) + 0(6U 2 ), 



one finds to first order 



.dC$(*i) fdu 



— i(uj+e a -eB)ti 

2tt 



A/^ 



Integrating this equation and using the boundary condition 

(i) 

*/3 



Cg(*->-oo) = 



one gets 



e -<(«+e a -efl)i 



(/3| ( 5f/H|a)+£rf QiV H(A| ( 5C/H| / i) 

A// 



(B3) 



27T oj + e a — Ep 

We finally arrive at the time-dependent Hartree-Fock single-particle wavefunctions up to the first order of 5U 

' doj e~ lut 



(B4) 



(B5) 



(B6) 



(B7) 



(B8) 



(B9) 



ct> a (rt;SU) = e- i£ " t \ ^ a {r) + £ 







2ir oj + e a — ep 



(/Wo;)|a> +J2T^(uj){X\SU(uj)\f J ,) 



doj e 1 ' 



2ir oj + e a — ep 



(P\6U(u)\a)* +£r£; A >)<A|5[/H| / z)' 



ipp(r) 

*m}- (bio) 



The second terms of these equations are respectively the va riations 8^j a (f,t) and J-0*(f, t). 
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